In [ ]:
# I load the needed libraries
library(dplyr)
library(scales)
library(GoFKernel)

library(mvtnorm)
library(gplots)

options(warn=-1)
Attaching package: ‘dplyr’


The following objects are masked from ‘package:stats’:

    filter, lag


The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


Loading required package: KernSmooth

KernSmooth 2.23 loaded
Copyright M. P. Wand 1997-2009


Attaching package: ‘gplots’


The following object is masked from ‘package:stats’:

    lowess


Preparation of the simulation¶

I load the functions from the class file:

In [ ]:
source("class_MCMC.R")

I define then the function that I want to use as output of the MCMCs:

In [ ]:
# Function to sampled from: n-dim gaussian with chosen sigmas and centers
posterior_g_inhom = function (theta) {

    sigmas = c(1:length(theta))
    centers = c(seq(length(theta), 1))

    product = 1
    for (i in 1:length(theta)) {
        product = product * exp(-(theta[i] - centers[i])**2/sigmas[i]**2)
    }

    return (product)

}

chosen_function = posterior_g_inhom

Then I only have to determine the parameters for the initialization = the "hyperparameters" of the simulations

In [ ]:
# The initial parameters are:
init = c(1, 1, 1)
std = diag(0.01, 3)

N = as.integer(1e5)
burn_in = as.integer(1e4)
# N = as.integer(1e4)
# burn_in = as.integer(1e3)

print_step = as.integer(1e2)
# print_init = as.integer(1e3)

N_tot = N + burn_in

# For Haario:
epsilon = 0.001

Simulations¶

In [ ]:
# MVTNORM 

# Evaluate then the MCMC
mcmc_g = random_steps_mvtnorm (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE)

# Selecting the sequence after the burn-in
mcmc_g = mcmc_g[burn_in:N, ]

# Plotting the results
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  94 %
In [ ]:
# MVTNORM GIBBS

mcmc_g = random_steps_mvtnorm_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  96.91909 %
In [ ]:
# # SIMPLE ADAPTIVE

# mcmc_g = random_steps_simple (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
#                                 gamma_function = gamma_series_exp, halved_step = burn_in)

# mcmc_g = mcmc_g[burn_in:N, ]

# show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
In [ ]:
# # SIMPLE ADAPTIVE GIBBS

# mcmc_g = random_steps_simple_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
#                                 gamma_function = gamma_series_exp, halved_step = burn_in)

# mcmc_g = mcmc_g[burn_in:N, ]

# show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
In [ ]:
# HAARIO

mcmc_g = random_steps_haario (func_wanted = chosen_function, theta_init = init, n_samples = N_tot,
                                sigma = std, print_accept=TRUE, t_0 = burn_in, eps = epsilon)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  16.84 %
Final mean =  2.994646 2.047772 0.9126068 
Final covariance matrix = 
          [,1]      [,2]      [,3]
[1,] 31.576237 21.650030  7.056395
[2,] 21.650030 18.945465  4.562009
[3,]  7.056395  4.562009 10.810059
In [ ]:
# HAARIO GIBBS

mcmc_g = random_steps_haario_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot,
                                    sigma = std, print_accept=TRUE, t_0 = burn_in, eps = epsilon)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  20.45485 %
Final mean =  2.988557 1.989517 0.8476119 
Final covariance matrix = 
          [,1]      [,2]     [,3]
[1,] 31.278268 19.636973 6.395148
[2,] 19.636973 15.870224 4.294283
[3,]  6.395148  4.294283 8.958365
In [ ]:
# RAO

mcmc_g = random_steps_AM_rao (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
                                gamma_function = gamma_series_exp, halved_step = burn_in/2)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  53.49 %
Final mean =  3.006974 2.036448 0.9892942 
Final covariance matrix = 
             [,1]         [,2]         [,3]
[1,]  0.512740395 -0.004800515 -0.005604978
[2,] -0.004800515  2.040432407  0.076991233
[3,] -0.005604978  0.076991233  4.499651563
In [ ]:
# RAO GIBBS

mcmc_g = random_steps_AM_rao_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
                                gamma_function = gamma_series_exp, halved_step = burn_in/2)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  64.85424 %
Final mean =  3.006909 2.013894 1.076841 
Final covariance matrix = 
             [,1]        [,2]        [,3]
[1,]  0.591654757 0.001114013 -0.02806483
[2,]  0.001114013 1.888830556  0.04854427
[3,] -0.028064831 0.048544275  3.87236218
In [ ]:
# GLOBAL

mcmc_g = random_steps_global (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
                                gamma_function = gamma_series_exp, halved_step = burn_in)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  99.15455 %
Final mean =  2.978676 0.8515176 1.301308 
Final lambda =  -1.182391 
Final covariance matrix = 
               [,1]           [,2]           [,3]
[1,]  4.940656e-324  4.940656e-324 -4.940656e-324
[2,]  4.940656e-324  4.940656e-324 -4.940656e-324
[3,] -4.940656e-324 -4.940656e-324  4.940656e-324
In [ ]:
# GLOBAL GIBBS

mcmc_g = random_steps_global_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
                                gamma_function = gamma_series_exp, halved_step = burn_in)

mcmc_g = mcmc_g[burn_in:N, ]

show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate =  94.61818 %
Final mean =  3.026764 1.83856 1.422002 
Final lambda =  -3.193156 
Final covariance matrix = 
           [,1]       [,2]       [,3]
[1,]  0.3821043 -0.0152492 -0.2344517
[2,] -0.0152492  2.2027019 -0.1621688
[3,] -0.2344517 -0.1621688  3.6836446